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Abstract 

A systematic comparison is conducted for pairing properties of finite systems at nonzero tem- 
perature as predicted by the exact solutions of the pairing problem embedded in three principal 
statistical ensembles, as well as the unprojected (FTBCSl+SCQRPA) and Lipkin-Nogami pro- 
jected (FTLNl+SCQRPA) theories that include the quasiparticle number fluctuation and coupling 
to pair vibrations within the self-consistent quasiparticle random-phase approximation. The nu- 
merical calculations are performed for the pairing gap, total energy, heat capacity, entropy, and 
microcanonical temperature within the doubly-folded equidistant multilevel pairing model. The 
FTLNl+SCQRPA predictions agree best with the exact grand-canonical results. In general, all 
approaches clearly show that the superfluid-normal phase transition is smoothed out in finite sys- 
tems. A novel formula is suggested for extracting the empirical pairing gap in reasonable agreement 
with the exact canonical results. 
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I. INTRODUCTION 



Pairing correlations are a fundamental feature responsible for the superconducting (su- 
perfluid) properties in many-body systems ranging from very large ones as neutron stars 
to tiny ones as atomic nuclei. In macroscopic systems such as superconductors, pairing 
correlations are destroyed as the temperature T increases and completely vanish at a value 
Tc ~ 0.57A(0), which is the critical temperature of the phase transition from the superfluid 
state to the normal one. Here A(0) is the value of the pairing gap at T = 0. The recent years 
witness a renewed interest in pairing correlations, which is supported by the exact solutions 
of the pairing problem in practice, the studies of unstable nuclei, the BCS to Bose-Einstein 
condensation crossover, etc. 

The superfluid properties and superfiuid-normal (SN) phase transition of infinite systems 
are accurately described by the Bardeen-Cooper-Schrieffer theory [l], where the average 
within the grand canonical ensemble (GCE) is used to obtain the occupation number in 
the form of Fermi-Dirac distribution for noninteracting fermions. The GCE consists of 
identically prepared systems in thermal equilibrium, each of which shares its energy and 
number of particles with an external heat bath. As compared to the other two principal 
thermodynamic ensembles, namely the canonical ensemble (CE) and the microcanonical one 
(MCE), the GCE is, perhaps, the most popular in theoretical studies of systems at finite 
temperature because it is very convenient in the calculations of average thermodynamic 
quantities ^. The CE is also in contact with the heat bath, but the particle number is the 
same for all systems. The MCE is an ensemble of thermally isolated systems sharing the 
same energy and particle number. In thermodynamics limit (i.e. when the system's particle 
number N and volume V approach infinity, but N/V is finite), fluctuations of energy and 
particle number are zero, therefore three types of ensembles offer the same average values for 
thermodynamic quantities. Thermodynamics limit works quite well in large systems as well 
where these fluctuations are negligible. The discrepancies between the predictions by three 
types of ensembles arise when thermodynamics is applied to small systems such as atomic 
nuclei or nanometer-size clusters. These systems have a fixed and not very large number of 

^ For instance, the double-time Green's functions G{t,t'), which are defined by using the GCE depend 
only on the time difference (t — t'), which greatly simplify the derivations in many statistical physics 
applications. 
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particles, their single-particle energy spectra are discrete with the level spacing comparable 
to the pairing gap. Under this circumstance, the justification of using the GCE for these 
systems becomes questionable. A number of theoretical studies have also shown that, in 
these tiny sy stems j_t her mal fiuctuations become so large that they smooth out the sharp SN 

1, y, 1^. As the result, the pairing gap never collapses but decreases 



transition 13, 14,151,161, 
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monotonously with increasing T. These predictions are in qualitative agreement with the 
results obtained by averaging the pairing energy in the CE built on the eigenvalues obtained 
by exactly solving the pairing problem 111, ll2|. As a matter of fact, even at T = 0, the 



exact pairing solution in nuclei shows a sizable pairing energy in the region where the BCS 
solution collapses [12]. In the literature so far, under the pretext that a nucleus is a system 
with a fixed number of particle, the thermodynamic averages of the exact solutions of the 
pairing Hamiltonian are usually carried out within the CE, and the results are compared 
with those obtained within different theoretical approximations at T 7^ 0. The latter, such 
as the BCS, Hartree-Fock-Bogoliubov (HFB) theories, etc., as a rule, are always derived 
within the GCE, where both energy and particle number fiuctuate. On the other hand, the 
well-known argument that the nuclear temperature should be extracted from the MCE of 
thermally isolated nuclei is also quite often debated and studied in detail j?]. 

These results suggest that a thorough comparison of the predictions offered by the exact 
pairing solutions averaged within three principal thermodynamic ensembles, and those given 
by the recent microscopic approaches, which include fluctuations around the thermal pairing 
mean field in nuclei, might be timely and useful. This question is not new, but the answers 
to it have been so far only partial. Already in the sixties Kubo pLSi] drew attention to the 
thermodynamic effects in very small metal particles. Later, Denton et al. [l4| used Kubo's 
assumption to study the difference between the predictions offered by the GCE and CE 
for the heat capacity and spin susceptibility within a spinless equidistant level model for 
electrons. Very recently, the predictions for thermodynamic quantities such as total energy, 
heat capacity, entropy, and microcanonical temperature within three principal ensembles 
were studied and compared in Ref. Iisj] by using the exact solutions of an equidistant 
multilevel model with constant pairing interaction parameter. However, no results for the 
pairing gaps as functions of temperature are reported. 

In the present paper, we carry out a systematic comparison of predictions for nuclear pair- 
ing properties obtained by averaging the exact solutions in three principal thermodynamic 
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ensembles as well as those offered by recent microscopic approaches to thermal pairing. 
For the latter, we choose the unprojected and particle-number projected versions of the 
FTBCSl+SCQRPA, which we recently developed in Ref. j^. This approach has been rigor- 
ously derived based on the same variational procedure used for the derivation of the standard 
BCS theory, taking into account the effects due to quasiparticle-number fluctuation (QNF) 
and coupling to the self-consistent quasiparticle random-phase approximation (SCQRPA) 
in addition to those of the QNF. 

The paper is organized as follows. The pairing Hamiltonian, its diagonalization, en- 
semble treatments of the exact pairing solutions, the main results of the FTBCSl and 
FTBCSH-SCQRPA theories, as well as their particle-number projected versions, are sum- 
marized in Sec. II. The numerical results obtained within the doubly-folded equidistant 
multilevel model llj are analyzed in Sec. III. The paper is summarized in the last section. 



where conclusions are drawn. 



II. EXACT SOLUTION OF PAIRING HAMILTONIAN 
A. Exact solution at zero temperature 

We consider a system of particles with single-particle energies ej, which are generated 
by particle creation operators on jth orbitals with shell degeneracies 2Qj {Qj = j + 1/2), 
and interacting via a monopole-pairing force with a constant parameter G. This system is 
described by the well-known pairing Hamiltonian 

H = J2'^N,-Gj2Ppf ^ (1) 

jm jj' 

where the particle-number operator Nj and pairing operator Pj are given as 

= J2 4-%- ' PhY. ^™«]- ' P^ = (^^)^ ' (2) 

m m>0 

with the symbol ~ denoting the time-reversal operator, namely ajf^ = {—)^~"^aj-m- For a 
two-component system with Z protons and neutrons, the sums in Eq. ([1]) run also over 

all jr^r, ir'^^'r-i ^"^^ wlth T = (Z, A^) . 

The pairing Hamiltonian ([T]) was solved exactly for the first time in the sixties by Richard- 



son 11]. By noticing that the operators J| = {Nj — flj)/2, = Pj , and Jj = Pj close 



an SU(2) algebra of angular momentum, the authors of Ref. [12] have reduced the problem 
of solving the Hamiltonian ([1]) to its exact diagonalization in the subsets of representations, 
each of which is given by a set basis states \k) = {Nj}) characterized by the partial 

occupation number Nj = (JJ + Qj)/2 and partial seniority (the number of unpaired parti- 
cles) Sj = Qj — 2Jj of the jth single-particle orbital. Here Jj{Jj + 1) is the eigenvalue of 
the total angular momentum operator {JjY ^ ^/"^j" + '^ji'^j ~ '^^^ partial occupation 
number Nj and seniority Sj are bound by the constraints of angular momentum algebra 
Sj < Nj < 2Qj — Sj with Sj < Qj. In the present paper, we use this exact diagonalization 
method to find the eigenvalues £s and eigenstates (eigenvectors) \s) of Hamiltonian ([T]). 
Each sth eigenstate \s) is fragmented over the basis states \k) according to |s) = Cjf'^\k) 
with the total seniority s = J2j Sj, and degenerated by 

(2^)! (2^)! 



4=n 



(3) 



_Sj\{2n - Sj)\ {sj - 2)!(2fi - Sj + 2)! 
where (C^'^'')^ determine the weights of the eigenvector components. The state-dependent 

(s) 

exact occupation number /j on the jth single-particle orbital is then calculated as the 
average value of partial occupation numbers Nj''^ weighted over the basis states \k) as 



B. Exact solution embedded in thermodynamic ensembles 

The properties of the nucleus as a system of interacting fermions at energy S can be 
extracted from its level density p{S, N) 

p(£,iV) = 5^<5(^-^i"))5(iV-n) , (5) 

where 8s^^ are the energies of the quantum states \s) of the n-particle system. Applying 
to the pairing problem, these energies are the eigenvalues, which are obtained by exactly 
diagonalizing the pairing Hamiltonian ([T]), as has been discussed in the previous section. 



1. Grand canonical and canonical ensembles 



Since each system in a GCE is exchanging its energy and particle number with the heat 
bath at a given temperature T = both of its energy and particle number are allowed 
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to fluctuate. Therefore, instead of the level density ([5]), it is convenient to use the average 
value, which is obtained by integrating ([5]) over the intervals in £ and A^. This Laplace 



transform of the level density p{E, N) deflnes the grand partition function A) |16| . 

Z{(3, A) = // p{S, N)e-^^^-^^^dNd£ = V e^^"Z(/3, n) , (6) 

where Z{l3, n) denotes the partition function of the CE at temperature T and particle number 
fixed at n, namely 

Z{I3, n) = Y,PsiP, n) , p,{(3, n) = Sr^e-^''"' . (7) 

s 

Within the GCE, the chemical potential A should be chosen as a function of T so that the 
average particle number (A^) of the system always remains equal to A^. The summations 
over n and s in Eqs. ([6]) and (I7j) are obtained by using Eq. ([5]) taking into account the 
degeneracy rfi"'' ([3]) of each sth state in the n-particle system, and carrying out the double 
integration with S functions. 

The thermodynamic quantities such as total energy S, heat capacity C within the GCE 
and CE are as usual given as 

where a=GC for the GCE and a=C for the CE. 

The thermodynamic entropy S^^^ is calculated within the GCE or CE based on the general 
definition of the change of entropy (by Clausius): 

dS = (3dS . (9) 

By using the differentials of lnZ(/?)a, one obtains 

4°^) = (3{{S)gc - AAT) + lnZ(/3, A) , = (3{S)c + lnZ(/5, N) . (10) 

The occupation number fj on the jth single-particle orbital is obtained as the ensemble 

(s) 

average of the state-dependent occupation numbers /j , namely 

,(GC) _ 1 V- As,n),n -p{ei"^~Xn) AC) _ ^ V- f^^'^^rj^p-^S^ 1 1 

for the GCE and CE, respectively. 



2. Microcanonical ensemble 



Different from the GCE and CE, there is no heat bath within the MCE. A MCE consists 
of thermodynamically isolated systems, each of which may be in a different microstate 
(microscopic quantum state) \s), but has the same total energy S and particle number N. 
Since the energy and particle number of the system are fixed, one should use the level density 
([5]) to calculate directly the entropy by Boltzmann's definition, namely 

S^^^\S) = \nn{S) , (12) 

where Q{S) = p{S)AS is the statistical weight, i.e. the number of eigenstates of Hamiltonian 
([1]) within a fixed energy interval AS. The condition of thermal equilibrium then leads to 



the standard definition of temperature within the MCE as 



16| 



^" dS ~ p{£) d£ ' ^ 

using which one can build a "thermometer" for each value of the excitation energy £ of the 
system. 

In practical calculations, to handle the numerical derivative at the right-hand side of Eq. 
f|T3|) . one needs a continuous energy dependence of p{£) in the form of a distribution in 
Eq. ([5]). This is realized by replacing the Dirac-(5 function 5{x) at the right-hand side of 
Eq. ([5]), where x = £ — £s, with a nascent 5-function 5a{x) (a > 0), i.e. a function that 
becomes the original 5{x) in the limit cr ^ ^. Among the popular nascent 5 functions 



are the Gaussian (or normal) distribution, Breit-Wigner distribution 21|, and Lorentz (or 
relativistic Breit-Wigner) distribution, which are given as 

5<t{.x)g = — ^e"i^ , 5a(x)Bw 



aV2TT TT + (T^ 



1 g£'^ 

respectively. In these distributions, cr is a parameter, which defines the width of the peak 
centered aX £ = £s- The full widths at the half maximum (FWHM) F of these distributions 



are F = 2crv2in2 ~ 2.36(j for the Gaussian distribution, and 2a for the Breit-Wigner and 



This replacement is equivalent to the folding procedure for the average density, discussed in Sec. 2.9.3 of 
the textbook taken at zero degree {M — 0) of the Laguerre polynomial l)1^{x). 
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Lorentz ones. The disadvantage of using such smoothing is that the temperature extracted 
from Eq. f|T3|) . of course, depends on the chosen distribution as well as the value of the 
parameter a. It is worth mentioning that changing a is not equivalent to changing AE" for 
the discrete p{£) used in calculating the statistical weight VL{£) in Eq. f[T2]) . since the wings 
of any distributions in Eq. (IT^ extend with increasing cr, whereas for the discrete spectrum, 
no more levels might be found in the low £ by enlarging A£. 

The generalized form of Boltzmann's entropy f|T2l) is the definition by von Neumann 

S = -Tr[plnp] , (15) 

which is the quantum mechanical correspondent of the Shannon's entropy in classical theory. 
By expressing the level density p{£) in Eq. ([5]) in terms of the local densities of states 

k s 

the entropy becomes 

k 

By using Eqs. ( |T3l) and ( fT7|) one can extract a quantity Tg as the "microcanonical temper- 
ature of each eigenstate s" . 



C. Determination of pairing gaps from total energies 



1. Ensemble-averaged pairing gaps 

Although the exact solution of Hamiltonian ([1]) does not produce a pairing gap per se, 
which is a quantity determined within the mean field, it is useful to define an ensemble- 
averaged pairing gap to be closely compared with the gaps predicted by the approximations 
within and beyond the mean field. In the present paper, we define this ensemble-averaged 
gap Aq, from the pairing energy £^'^1^ of the system as follows 

A.=v'^G^, = (^:)« - (^)i°)=2E^•[9-f/.^"^/f^ (18) 

within the GCE (a = GC), CE (a = C), and MCE (a = MC), where for the latter we 
put £^p"j^j. = £paiT:{s) with (£^)mc = £s^^ ■ The term {£)a^ denotes the contribution from the 
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energy 2 Qjejfj"'^ of the single-particle motion described by the first term at the right- 
hand side of Hamiltonian ([1]), and the energy —G of uncorrelated single-particle 
configurations caused by the pairing interaction in Hamiltonian ([1]). Therefore, subtracting 
the term ^ from the total energy {S)a yields the residual that corresponds to the energy 
due to pure pairing correlations. By replacing fj""^ with f|, one recovers from Eq. ( fT8|) the 
expression S^^^^ = — Al^^g/G of the BCS theory. Given several definitions of the ensemble- 
averaged gap existing in the literature, it is worth mentioning that the definition (ITSl) is 
very similar to that given by Eq. (52) of Ref. [l8|, whereas, even within the CE, the gap 



Ac is different from the canonical gap defined in Refs. 19|, l20|], since in the latter, the term 



{S)q^ is taken at G = 0. The pairing energy S^'^jj. in Eq. (fT8|) is also different from the 
simple average value ~G^--,{P'-Pji)a of the last term of Hamiltonian ([T]) as the latter still 
contains the uncorrelated term — G 



2. Empirical determination of pairing gap at finite temperature 

The simplest way to empirically determine the pairing gap of a system with particles 
in the ground state (at T = 0) is to use the so-called three-point formula A'^^)(A^), which is 
given by the odd-even mass difference between the ground-state energies of the A^-particle 



system and the neighboring systems with ± 1 particles 
of this formula to T 7^ reads 



A straightforward extension 



A(3)(;3, N) ^ ^-^[{SiN + 1)), - 2{£iN))^ + {SiN - 1)) J , (19) 



which was used, e.g., in Ref. 2^ to extract the thermal pairing gaps in wolfram and 
molybdenum isotopes. The four-point formula represents the arithmetic average of the 
three-point gaps over in the neighboring systems with A^ and A^ — 1 particles, namely 

A(^) (/3, N) = ^ [ A(3) (A^) + A(3) (A^ - 1)] . (20) 

A drawback of the gaps A^^\P,N) and A^^\P,N), defined in this way, is that they still 
contain the admixture with the contribution from uncorrelated single-particle configurations. 
The later increases with increasing T. Therefore, Eq. (fT9|) and, consequently, Eq. (120|) do 
not hold at finite temperature. To remove this contribution so that the experimentally 
extracted pairing gap is comparable with Aq, in Eq. (ITSl) . we propose in the present paper 



an improved odd-even mass difference formula at T 7^ as follows. Using Eq. (fTSll to express 
the total energy {S{N))a of the system in terms of Aa{N) and (£^(A^))°, we obtain 

{sm. = i^m^^^ - ^^§^ ■ (21) 

where {S{N))a is the experimentally known total energy of the system with particles 
at T 7^ 0, whereas A{l3,N) is the pairing gap of this system to be determined. Replacing 
{S{N))a in the definition of the odd-even mass difference (fT9l) with the right-hand side of 
Eq. ( l2Ti) . we obtain a quadratic equation for the three-point A^^\P, N) 

A(3)(/?, N) = {-ifi^lmN + 1))„ + {S{N - 1)U - {S)l + . (22) 

The discriminant of this equation is equal to G\Jl — AS'/G, where 

S' = l [{SiN + 1)U + {S{N - 1)) J - (^(iV)) . (23) 

Therefore the condition for Eq. (l22l) to have real solutions is S' < G/A. Including both 
cases with even and odd A^, the positive solution of Eq. fl2^ is always possible provided 
5" < 0, which reads 

n[ I ^1 

(24) 



-i)^ + a/i-4 



The quantity 5" differs from the conventional odd-even mass difference shown in the square 
brackets of f|T9l) by the contribution due to the uncorrelated single-particle motion, i.e. the 
last sum containing G in the definition of {S)a^ in Eq. (ITSl) . The latter is zero only at G = 
0, which yields A^^\P,N) = 0, as can be seen from Eq. fl24l) as well. In this case both S' 
and the expression in the square brackets of Eq. (1191) vanish as they are just the difference 
of the Hartree-Fock energies {E{N + l))i°^ + {E{N - l))i°^ - 2{E{N))l^\ which is zero at 
G = 0. Moreover, while the odd-even mass difference of Eq. f|T9|) can be positive or negative 
depending on whether is even or odd, S' should be always negative as discussed above. 
The gap A^^\P,N) extracted from Eq. fl24l) is, therefore, consistent with the result of the 
exact calculation at zero temperature, where the pairing gap is zero only at G = 0, and 



increases with G (See e.g. Fig. 1 - (a) of Ref. [15|]). As compared to the simple finite- 
temperature extension of the odd-even mass (1191) . the modified gap A^^\P,N) is closer to 
the ensemble-averaged gap Aq,(A^) (fT8|) since it is free from the contribution of uncorrelated 
single-particle configurations. In Eq. (|24|) . the energies {£{N + 1))^ and {S{N — l))a can be 
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extracted from experiments, whereas the pairing interaction parameter G can be obtained 
by fitting the experimental values of A(T = 0, A^). The energy {£{N))^a'^ remains the only 
mo del- dependent quantity being determined in terms of the single-particle energies and 
single-particle occupation numbers /j"^- As a matter of fact, since the pairing gap A*^'^^ {(3, N) 
of the A^-particle system at the left-hand side of the expression (IT^ is also present in the 
total energy {S{N))a of the same system at the right-hand side of (fTOl) . the former is simply 
extracted from the latter by using Eqs. f|T8l) and fl2Tl) . As the result, the modified gap 
A^^\P,N) of the system with A^ particles explicitly depends now on {SiN))'^^ of the same 
system rather than on its total energy {£(N)). The modified four-point gap A^^\l3,N) is 
then obtained from the modified three-point gaps A^^'>(P,N) and A^^\P,N — 1) by using 
the definition 



III. FTBCSl-FSCQRPA AND FTLNl-^SCQRPA 

The FTBCSl+SCQRPA includes the effects due to quasiparticle-number fiuctuation and 
coupling to the SCQRPA vibrations, which are neglected within the standard BCS theory. 
The derivation of the FTBCSl+SCQRPA has already been given and discussed in detail in 
Ref Q]. Therefore we give below only the main results, which are necessary to follow the 
numerical calculations in the present paper. 

The rigorous derivation of the FTBCSl+SCQRPA theory follows the standard variational 
procedure, which is used to derived the BCS theory. By using the Bogoliubov's canonical 



transformation 



23| 



the pairing Hamiltonian ([1]) is expressed in terms of quasiparticle operators, a]^ and a^m, as 
iJq.p., whose explicit form is given in many papers, e.g., Eqs. (3), (8) - (14) of Ref . (gj. The 
Uj and Vj coefficients of the Bogoliubov's transformation ( !25l) are determined by minimizing 
the GCE average value of the Hamiltonian Ti = -f^q.p. — XN. This leads to the equation 

{[H, A]])gc = , where A] = ^ • (^6) 

V^h j=i 

The final result yields the equation for the level-dependent pairing gap for a system with 
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even number of particles, 



Aj = A + SAj , A = G^Qf{l - 2nf)ufVf , SAj = 2G 



1 - 2n. 



with the QNF defined as 



and the equation for the average particle number N, 



n 



(27) 



(28) 



(29) 



The Uj and Vj coefficients in Eqs. ( 1271) and ( l29l) are determined as 



e' - Gv^. - A 



with the renormalized single-particle energies e'j. 

G 



Wiil- 2n, 



(30) 



(31) 



and the quasiparticle energies 



E, 



(32) 



For a system with an odd number of particles, the blocking effect caused by the unpaired 
particle should be taken into account in Eqs. (1271) and (12^ . In the present paper, for simplic- 
ity, we do not consider systems with odd particle numbers within the FTBCSl+SCQRPA 
and FTLNl+SCQRPA. 

The pair correlators {j^^j^-,^-) and {J^^Aji) in Eq. (131]) are determined by numerically 



solving a set of coupled equations (47) and (48) of Ref. [9j, which contain the and y'^ 
amplitudes of the SCQRPA equations. The details of the derivation of the SCQRPA are 



given in Ref. 



2J]. The SCQRPA equations are solved self-consistently with the gap and 



number equations and (12^ . The quasiparticle occupation numbers rij are then found 
by solving a set of equations that include coupling of quasiparticle density operators ajmaim 
to the SCQRPA phonon operators. The result yields the integral equation (69) of Ref. [9| 
for rij. To compare with the level-independent gap such as the BCS one, the level- weighted 

gap, 

A = Y,^,A,/Y,^, , (33) 
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is used instead of Eq. (1271) . The total energy (i^)FTBCsi+scQRPA is calculated by averaging 
the quasiparticle representation -f^q.p. of the Hamiltonian within the GCE, i.e. 

{£) = (//q.p.)GC , (34) 

and the heat capacity C is then found from Eq. ([8]). The thermodynamic entropy is obtained 
by integrating Eq. IQ 

5th = [ -Cdr. (35) 



The BCS equations [l|, l23| are recovered from Eqs. (1271) and (129|) by neglecting the QNF 
57V/ [i.e. 6Aj = in Eq. ([27])] together with the pair correlators {A^jA^'j/^j) and {A^jAj') in 
Eq. (1311) . and assuming nj to have the form of Fermi-Dirac distribution for noninteracting 
quasiparticles, i.e. setting 



1 



= = ef^^^ + 1 ■ ^^^^ 
The BCS is known to violate the particle number. This causes a certain quantal particle- 
number fluctuation around the average value determined by Eq. (l29l) even at T = 0. The FT- 
BCSl+SCQRPA takes into account only the thermal effect in terms of QNF but does 
not remove the quantal fluctuation, which is a feature inherent in the BCS wave functions. 
To cure this inconsistency, a proper particle-number projection (PNP) needs to be carried 



out. The Lipkin-Nogami (LN) method [251] is an approximated PNP before variation, which 
is widely used in nuclear study because of its simplicity. This method has been implemented 
into the FTBCSl+SCQRPA, and the ensuing approach is called the FTLNl+SCQRPA (See 
Sec. II. C. 2. of Ref. {9]). However, the LN method can approximately eliminate only the 
quantal fluctuations due to particle-number violation within the BCS theory. These quantal 
fluctuations are different from the thermal particle-number fluctuations, which always arise 
from the exchange of particles between the systems in the GCE. The LN method is, there- 
fore, not sufficient to remove the particle-number fluctuations within the GCE. To avoid the 
thermal particle-number fluctuations, the average in the CE should be used instead. Unfor- 
tunately, the methods of equilibrium statistical physics applied to the nuclear theories such 
as the Matsubara-Green's function and/or double-time Green's function techniques, which 
are used to derived the BCS and QRPA equations at finite temperature, are all based on 



the GCE. The comp 
projection operator 



ete particle-number projection based on the applying particle-number 
26] at finite temperature for these approaches still remains a subject 



under study [See, e.g. Ref. [2 
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In this paper, the numerical results obtained within both particle-number unprojected 
(FTBCSl+SCQRPA) and projected (FTLNl+SCQRPA) approaches will be compared with 
those given by averaging the exact pairing solution with the principal thermodynamic en- 
sembles. 

IV. NUMERICAL RESULTS 

A. Details of numerical calculations 

The schematic model employed for numerical calculations consists of particles, which 
are distributed over Q = N doubly-folded equidistant levels (i.e. with the level degeneracy 
2Qj = 2). These levels, whose energies are ej = e[j — {Q + l)/2] (j = 1, . . . , interact 
via the pairing force with a constant parameter G. The model is half-filled, namely, in 
the absence of the pairing interaction, all the lowest Q/2 levels (with negative single-particle 
energies) are filled up with A^ particles, leaving Q/2 upper levels (with positive single-particle 
energies) empty ^. It is worth mentioning that the extension of the exact solution to T 7^ 
is not possible at a large value of Q = N. For the present schematic model, the number of 
eigenstates ns, each of which is 2'^-degenerated, increases almost exponentially with Q 

^^(^) = CU + EC?xC;;;^_., (37) 
s 

where = m\/[n\{m — n)\] and Npair = A^/2 is the numbers of pairs distributed over Q 
single-particle levels. The sum in Eq. (l37j) runs over all the values of total seniorities S 
= 0, 2, . . . ,Q . Therefore, a.t Q = N = 16 there are 5196627 states, which corresponds 
to the order ~ 2.7x10^'^ for the square matrix to be diagonalized. This makes the finite- 
temperature extension of the exact pairing solution practically impossible for Q > 16 since 
all the eigenvalues must be included in the partition function. Therefore, in the present 
paper, we limit the calculations np to Q = N = 14, for which there are 73789 eigenstates. 
For the GCE average with respect to the system with A^ particles and Q = N levels, the 
sum over particle numbers n runs from rimin = 1 to Umax = 2^2 — 1 with the blocking effect 
caused by the odd particle properly taken into account. The calculations are carried out by 

This model is also called Richardson's model, picket-fence model, ladder model, multilevel pairing model, 
etc. in the literature. 
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using the level distance e = 1 MeV and the pairing interaction parameter G = 0.9 MeV. 
With these parameters, the values of the pairing gap obtained at T = are equal to around 
3, 3.5, and 4.5 MeV for = 8, 10, and 12 in qualitative agreement with the empirical 
systematic for realistic nuclei |l^. 



B. Results within GCE, CE and BCS-based approaches 

Shown in Fig. [1] are the pairing gaps, total energies, and heat capacity obtained for = 
8, 10, and 12 as functions of temperature T within the FTBCS, FTBCSl+SCQRPA, and 
FTLNl+SCQRPA, along with the corresponding results obtained by embedding the exact 
solutions (eigenvalues) of the Hamiltonian ([1]) in the CE and GCE. These latter results 
using the exact pairing solutions are referred to as CE and GCE results hereafter. For the 
FTBCSl+SCQRPA and FTLNl+SCQRPA the level-weighted gaps A ([33]) are plotted in 
Figs. [I](al) -[I](a3). 

It is seen from this figure that the GCE results are close to the CE ones for the gaps, 
obtained from Eq. (fT8|) . as well as for the total energies, obtained by using Eq. ([8]). 
Among three systems under consideration, the largest discrepancies between the CE and 
GCE results are seen in the lightest system (A^ = 8), for which the GCE gap is slightly 
lower than the CE one, and consequently, the GCE total energy is slightly higher than that 
obtained within the CE. As A^ increases, the high-T values of the GCE and CE gaps become 
closer, so do the corresponding total energies. Different from the BCS results (dotted lines), 
which show a collapse of the gap and a spike in the temperature dependence of the heat 
capacity at T = Tc, no singularity occurs in the GCE and CE results. Both GCE and 
CE gaps decrease monotonously with increasing T and remain finite even at T ^ 5 MeV. 
The FTBCSl+SCQRPA and FTLNl+SCQRPA predictions for the pairing gap are found 
in qualitative agreement with the GCE and CE results [Figs. [1] (al) -[1] (a3)]. Because 
of a different definition of the pairing gap in the exact solutions within the GCE and/or 
CE, where actually no mean-field gap exists (See Sec. II CI), one cannot expect a more 
quantitative agreement between the predictions by the FTLNl+SCQRPA and the GCE 
(CE) results. In this respect, the modified gap (|24l) yields a better quantitative agreement 
with the GCE (CE) results, as will be seen later in Sec. IIVDI The two BCS-based approaches 
differ noticeably only at T < Tc, where the FTLNl+SCQRPA gap, due to PNP, practically 
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FIG. 1: (Color online) Pairing gaps A, total energies (if), heat capacities C and thermodynamic 
entropies, obtained for = 8, 10, and 12 {G = 0.9 MeV) within the FTBCS (dotted lines), 
FTBCSl+SCQRPA (thin dashed hnes), FTLNl+SCQRPA (thick dashed lines), CE (dash-dotted 
lines), and GCE (solid lines) vs temperature T. 

coincides with the GCE and CE results at T < 0.5 - 1 MeV. At T > the predictions 
by two approaches start to converge to the same value, which decreases with increasing T, 
remaining smaller than the GCE and CE gaps. 

For the same reason, at T < Tc, where the total energy predicted by the 
FTLNl+SCQRPA agrees very well with the GCE and CE results, the FTBCSl+SCQRPA 
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energy is significantly larger [Figs[T](bl) -[I](b3)]. At T > Tc, one finds a remarkable agree- 
ment between the energies predicted by the FTBCSl+SCQRPA, FTLNl+SCQRPA, and 
that obtained within the GCE. This seems to be a natural consequence, given the fact that 
the two BCS-based approaches are derived by using the variational procedure within the 
GCE. The energies do not suffer either from the difference in the definitions as the pairing 
gaps do, as has been discussed above. For the heat capacities [Figs[T](cl) -[I](c3)], the spike 
obtained at T = Tc within the FTBCS theory is completely smeared out within the GCE, 
CE as well as the FTBCSl+SCQRPA and FTLNl+SCQRPA, where only a broad bump is 
seen in a large temperature region between and 3 MeV. At T < 1 - 1.2 MeV, the difference 
between the GCE and CE energies leads to a significant discrepancy between the GCE and 
CE values for the heat capacity. As the FTBCSl+SCQRPA and FTLNl+SCQRPA heat 
capacities are close to the GCE values, this explains the discrepancy reported in Figs. 4 (b) 



and 4 (c) of Ref. [9] between these results and the predictions obtained within the CE by 



the quantum Monte-Carlo calculations, which use a model Hamiltonian with same monopole 
pairing interaction. The thermodynamic entropies Sth are shown as functions of T in Figs. 
[H (dl) -[1] (d3). Once again the FTLNl+SCQRPA results for the thermodynamic entropy 
Sth agree very well with obtained within the GCE, whereas such good agreement is seen 
for the FTBCS 1-SCQRPA results only at T < 1 MeV. The CE thermodynamic entropy 
is significantly lower than the values obtained in all other approaches under consideration, 
which are based on averaging within the GCE. On the other hand, the BCS theory strongly 
overestimates the thermodynamic entropy at T > Tc. 

C. Results within MCE 

The values of temperature within the MCE as extracted by using Eq. (|T3l) are plotted in 
Fig. [2] along with the CE results against the excitation energy S*. The CE definition of the 
latter is = {S(T))c — {S(T = 0))c- For comparison, the results obtained by using the 
entropy (ITTl) are also presented in the top panels [Figs. [2] (a) -[2] (c)]. They show the values 
of the eigenstate temperatures Tg, which scatter widely around the heat bath temperature 
(i.e. the CE result) [Figs. [2] (a) -[2] (c)]. Many of these values are even negative. Since 
the eigenstate temperatures T^ are related with the spread of the exact eigenfunctions over 
the "unperturbed" basis states k with the weights [C^*'']^, they need not follow the trend 
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FIG. 2: (Color online) Temperature extracted from Eq. ()13p within the MCE (dots) vs excitation 
energy £* in comparison with the CE results (dash-dotted line) for N = 8, 10, and 12. The results 
in the top panels, (a) - (c), are obtained from the definition of Neumann for entropy Eq. (jlSp . 
whereas those in the middle panels, (d) - (f), and bottom ones, (g) - (i), are calculated from 
Boltzmann's entropy Eq. (I12p with two different values of the energy interval AE in the statistical 
weight 

of the heat bath (or canonical) temperature, which depends just on the level density [Eqs. 
( |T2l) and ( |T3i) ]. In fact, with increasing the energy interval AE, within which the levels are 
counted, the values of the MCE temperature determined by Eq.f ll3l) gradually converge to 
the CE values [Figs. [2] (a) -[2](i)]. This means that, the thermal equilibrium within the MCE 
for the present isolated pairing model can be reached only at large and dense spectrum 
(small level spacing). Moreover, a full thermalization in a system with pure pairing is a 
subject under question. In Ref. 0, IsS], the thermalization of the system is characterized 
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FIG. 3: (Color online) Temperatures [(a) - (c)], entropies [(d) - (f)], and pairing gaps [(g) - (i)] 
within the MCE as functions of excitation energy E* for = 8 (G = 0.9 MeV) obtained by 
using the Gaussian, Lorentz, and Breit-Wigner distributions from Eq. (jl4p for the level density at 
different values of the parameter o. The dash-dotted lines show the CE values. 

by the single-particle temperature, which is obtained by fitting the occupation numbers of 
the individual eigenstates to those given by the Fermi-Dirac distribution. The numerical 
results Fig. 12 of Ref. 3] show that, the temperature extracted from the density of states 
in ^^^Sn by using Eq. (fT3l) agrees with the single-particle temperature only when all the 
residual interactions are taken into account. With the pure pairing interaction alone, the 
single-particle temperature shows a low temperature of the whole system, whereas the MCE 
temperature ( fT3|) is a hyperbola as a function of the excitation energy with a singularity in 
the middle of the spectrum, where it turns negative [See also Figs. 53 and 54 of Ref. [Tlj]. 

The values of MCE temperature, entropy and gap obtained by using the Gaussian, 
Lorentz, and Breit-Wigner distributions from Eq. f|T^ are shown in Fig. [3] as functions 
of excitation energy E* . While the fluctuating behavior of the microcanonical temperature 
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can be smoothed out by increasing the parameter a in all three distributions, we found that 
only the Gaussian distribution can simultaneously fit both the temperature and entropy 
[Figs. [3] (a) and [3] (d)]. The Lorentz distribution can fit only the MCE temperature to the 
CE one, but fails to do so for the entropy [Figs. [3] (b) and [3] (e)], whereas the Breit-Wigner 
distribution can fit the MCE temperature to the CE value only at high excitation energies 
[Figs. [3] (c)]. A similar result is seen for the pairing gaps as functions of £* ^ where the 
Gaussian fit gives the best performance among the three distributions [Compare Figs. [3] (g), 
|3] (h), and [3] (i)]^. We conclude that the Gaussian distribution should be chosen as the best 
one for smoothing the level density p{£) in Eq. (I5l) to extract the MCE temperature. 



D. Pairing gaps extracted from odd-even mass differences 

We extracted the pairing gaps A*^*)(/3, A^) (i = 3 and 4) by using the simple extension of 
the odd-even mass formula to T 7^ in Eq. ( |T9l) as well as the modified gaps A'^*)(/5,A^) 
from Eq. ^2M . and the canonical gaps Aq^ from Eq. ( fTSl) for several particle numbers up 
to = 12. In these calculations, the blocking effect caused by the unpaired particle in 
the systems with an odd particle number is properly taken into account in constructing the 
basis states when diagonalizing the pairing Hamiltonian. The results obtained for A^ = 9 
and 10 {Vt = 10) are displayed in Fig. HI First of all, the values for S' are found to be 
always negative, and increases with T > 1 MeV to reach a value of around —2 MeV at 
T = 5 MeV, and vanish at very high T. By comparing Figs. H] (a) and H] (b) one can see 



a clear manifestation of the parity effect 18|], which causes the reduction of the three-point 
gap the in the system with an odd particle number (A^ = 9) due to the blocking effect. For 
lights systems as those considered here, this reduction is rather strong (about 1 MeV at 
T = 0). With increasing T, as the thermodynamics weakens the effect of a single-unpaired 
particle, the parity effect starts to fade out in such a way that the three-point canonical 
gap {N = 9) slightly increases with T up to T ~ 1 MeV, starting from which it is even 
slightly larger than A^^\n = 9). This feature is found in qualitative agreement with the 



^ The non-vanishing values of T and S a,t £* — m Figs. [3] (a), [3] (c), [3] (d), [3] (f) are the artifacts due 
to the use of the Gaussian and Breit-Wigner distribution functions (14) to smooth out the discrete level 
density. These two distributions have non-zero values '--^ 1/cr at a; = 0, where the Lorentz distribution 
vanishes. 
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FIG. 4: (Color online) Pairing gaps extracted from the odd-even mass differences as functions of 
T for iV = 10 (a,c) and iV = 9 (b,d) {n = 10, G = 0.9 MeV). The thin solid and thick sohd lines 
denote the gaps A(*)(/?,iV) {i = 3, 4) from Eq. 1^, and the modified gaps AW(/?,iV) from Eq. 
()24p . respectively. The dash-dotted lines are the canonical results A(^\ The upper panels (a) and 
(b) show the three-points gaps {i = 3), whereas the the corresponding four-points gaps (i = 4) are 
shown in the lower panels (c) and (d). 



results obtained in Fig. 25 of Ref. [l&l for ultrasmall metallic grains. It is also seen in 
Fig. m that the naive extension of the odd-even mass formula to T ^ 0, resulting in the 
gap A(*)(/?, A^) (thin solid lines), fails to match the temperature-dependence of the canonical 
gap Aq^ (dash-dotted lines). The former even increases with T at T < 1 MeV, whereas 
it drastically drops at T > 1 - 1.5 MeV, resulting in a very depleted tail at T > 2 MeV 
as compared to the canonical gap A^^"*. Moreover, the three-point gap A*^'^^(/?, = 9) even 
turns negative at T > 2.4 MeV, suggesting that such simple extension of the odd-even mass 
difference to finite T is invalid. At the same time, the modified gap A^'^\P,N) (thick solid 
line) given by Eq. ( 12^ is found in much better agreement with the canonical one. At T < 
1.5 MeV, the three point-gap A*^^^(/3, A^) is almost the same as the canonical gap. At higher 
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T, it becomes larger (smaller) than the canonical value for the system with an even (odd) 
A^, however the systematic of the results of our calculations up to = 12 show we that this 
discrepancy decreases with increasing the particle number. The source of the discrepancy 
resides in the assumption of the odd-even mass formula that the gap obtained as the energy 
difference between the systems with + 1 and particles is the same as that obtained from 
the energy difference between systems with and A^ — 1 particles. This assumption does 
not hold for small A^ (Cf. Ref. 29|). The average in the four-point gap nearly eliminates 
this difference. As a result, the modified four-point gaps A^^\P, N) practically coincide with 
the canonical gaps. A natural consequence of the average in the definition of the four-point 
gap fl2Ul) is that the gaps obtained in the systems with A^ and A^ — 1 particles are now nearly 

nnnnn; 

the same. The pairing gaps predicted by a number of alternative theories p|, 15|, [g, 18|, 
including those discussed in the present work, are in closer agreement to the GCE and CE 
gaps rather than to the gap A{f3,N) from Eq. (fT9|) . Therefore, the comparison in Fig. H] 
suggests that formula (j24l) is a much better candidate for the experimental gap at T 7^ 0, 
rather than the simple odd-even mass difference (ITQll . 



V. CONCLUSIONS 



In the present work, a systematic comparison is conducted for pairing properties of finite 
systems at finite temperature as predicted by the exact solutions of the pairing problem 
embedded in three principal statistical ensembles, as well as by the recently developed 
FTBCSl (FTLN1)+SCQRPA. The analysis of numerical results obtained within the doubly- 
folded equidistant multilevel model for the pairing gap, total energy, heat capacity, entropies, 
and MCE temperature allows us to draw following conclusions. 

1) The sharp SN phase transition is indeed smoothed out in exact calculations within all 
three principal ensembles. The results obtained within GCE and CE are very close to each 
other even for systems with small number of particles. As for the MCE, although it can 
also be used to study the pairing properties of isolated systems at high-excitation energies, 
there is a certain ambiguity in the temperature extracted from the level density due to the 
discreteness of a small-size system. This ambiguity, therefore, depends on the shape and pa- 
rameter of the distribution employed to smooth the discrete level density. We found that, in 
this respect, the normal (Gaussian) distribution gives the best fit for both of the temperature 
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and entropy to the canonical values. The wide fluctuations of MCE temperature obtained 
here also indicate that thermal equilibrium within thermally isolated pure-pairing systems 
might not be reached. On the other hand, it opens an interesting perspective of studying the 
behavior of phase transitions in finite systems within microcanonical thermodynamics [s^ 
by using the exact solutions of pairing problem. 

2) The predictions by the FTBCSl+SCQRPA and FTLNl+SCQRPA are found in rea- 
sonable agreement with the results obtained by using the exact solutions embedded in the 
GCE and CE. The best agreement is see between the FTLNl+SCQRPA and the GCE 
results. Once again, this is a robust confirmation that quasiparticle-number fluctuation, 
included in these approximations, is indeed the microscopic origin of the strong thermal 
fluctuations that smooth out the sharp SN phase transition in finite systems. 

3) We suggest a novel formula to extract the pairing gap at finite temperature from the 
difference of total energies of even and odd systems where the contribution of uncorrelated 
single-particle motion is subtracted. The new formula predicts a pairing gap in much better 
agreement with the canonical gap than the simple finite-temperature extension of the odd- 
even mass formula. 
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